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In this work we extend the applicability of the microcanonical ensemble simulation method, 
originally proposed to study the Ising model (A. Hiiller and M. Pleimling, Int. Journal of Modern 
Physics C, 13, 947 (2002)), to the case of simple fluids. An algorithm is developed by measuring 
the transition rates probabilities between macroscopic states, that has as advantage with respect to 
conventional Monte Carlo NVT (MC-NVT) simulations that a continuous range of temperatures are 
covered in a single run. For a given density, this new algorithm provides the inverse temperature, 
that can be parametrized as a function of the internal energy, and the isochoric heat capacity is then 
evaluated through a numerical derivative. As an illustrative example we consider a fluid composed 
of particles interacting via a square-well (SW) pair potential of variable range. Equilibrium internal 
energies and isochoric heat capacities are obtained with very high accuracy compared with data 
obtained from MC-NVT simulations. These results are important in the context of the application 
of Hiiller-Pleimling method to discrete-potential systems, that are based on a generalization of the 
SW and Square-Shoulder fluids properties. 


PACS numbers: 05.10.-a, 05.20.Jj, 51.30.-H 


I. INTRODUCTION 


The random cluster model was first introduced by 
Kasteleyn and Fortuin 0,0 , and by Coninglio and Klein 
[aj for lattice systems. Swendsen and Wang 0 0 devel¬ 
oped the cluster algorithm based on the physics of this 
model and were able to implement an algorithm using 
the Widom-Rowlinson and Stillinger-Helfand models for 
fluid mixtures 00, an example that numerical simula¬ 
tion methods first developed for lattice systems can also 
be extended to off-lattice systems. 

Hiiller and Pleimling 0 have developed a success¬ 
ful cluster algorithm to study the Ising model in two 
and three dimensions within the microcanonical ensem¬ 
ble based upon the Broad Histogram Method (BHM) 0 - 
0. The algorithm works very well because in the Ising 
model the energy and magnetization have discretized val¬ 
ues. The Hiiller-Pleimling method (HPM) has two im¬ 
portant advantages compared with the method proposed 
by Creutz (CM) 0]: First, in a single run a wide range 
of energies can be covered, instead of a fixed one as in 
CM . Second, the heat capacity can be obtained with 
a numerical derivative, whereas for CM it is necessary 
to calculate the fluctuations in the Fourier transform of 
the energy. In the CM approach, however, no floating¬ 
point operations are required and no random numbers 
are needed at all if sequential updating is used. 
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In terms of efficiency the HPM approach compares very 
well with the method proposed by Wang and Landau 
(WLM) 0, that covers a wide range of temperatures 
too. For example, for the 32 x 32 Ising model, the average 
errors for the energy are 0.025% and 0.035% for HPM 
and WLM, respectively, considering the same number of 
Monte Carlo updates (7 x 10 5 , see reference 0 ). 

The first successful attempt to extend BMH to contin¬ 
uous systems was made by Munoz and Herrmann [0 . 
In this work we use HPM to study simple fluids, taking 
advantage that in this method the entropy is evaluated 
as a function of the energy. The extension presented 
here can be applied to discrete-potential (DP) systems, 
that also have a discrete set of energy values as in the 
Ising model. The SW potential is a very simple DP sys¬ 
tem after the hard-sphere potential, whose structural and 
thermodynamic quantities have been very well character¬ 
ized along the years fl5l - f27| . The SW system has been 
extensively applied within different statistical mechanics 
approaches 0, [20] , since is usually easier to implement 
than other pair potentials in statistical-mechanics theo¬ 
ries, as the Discrete Potential and Multipolar Discrete 
Perturbation Theories 0, 0] , as well as the Statistical 
associating fluid theory for chain molecules with attrac¬ 
tive potentials of variable range (SAFT-VR) 0]. These 
and other theories have made possible to describe phase 
diagrams of real complex fluids. 

The SW potential is a radial potential defined as 

{ oo if r < a 

—e if a < r < Act , (1) 

0 if r > Act 

where r is the separation distance between the centers of 
two particles; er represents the particle’s hard-core diam- 
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eter; A is the potential range of an attractive interaction 
of depth -e. 

In section II the method proposed in this work is 
described in detail, and its application to SW systems 
is presented in Section III. Among the thermodynamic 
quantities that can be calculated, we focus our attention 
on the internal energy and the isochoric heat capacity, 
since these quantities are more sensitive when perturba¬ 
tion theories are used. As we shall see in this section, 
both properties are predicted very accurately with mi- 
crocanonical simulations when compared with MC-NVT 
data. Finally, in Section IV the main conclusions of this 
work are given. 


II. SIMULATION METHOD 


In order to describe the methodology we consider a 
system of N particles confined within a volume V. The 
available energy levels will be labeled as A 1 ,, defined as 


Ei = -ie= ^2 ( 2 ) 

k,l^k 


where r^i is the distance between the centers of the k 
and l particles and i is the number of pair of particles 
that satisfies a < r^i < Act. We denote by Q,(E{) the 
number of configurations, or microstates, that share the 
same energy Ei. 

A new microstate is generated when a given mechanism 
is applied to a previous microstate; for a fluid composed 
of spherical particles the mechanism that has to be con¬ 
sidered is a displacement of a particle, Ar. When Ar 
is applied to all N particles for each tt(Ei) microstates, 
then Nfl(Ei) new microstates are created, but only a 
small number will have energy Ej, that will be denoted 
as Vij ; in order to satisfy the reversibility condition, it 
is required that Vij =Vji, a relation that holds for both 
continuous and discontinuous systems, since the BHM 
method is an exact theory (see reference [H|). The set of 
movements counted applying BHM is different to those 
followed to construct the sequence of visited, averaged 
states, that guarantee equiprobability. Consequently, it 
is necessary to consider two different protocols of allowed 
movements, depending on the reversibility or equiproba¬ 
bility conditions to be satisfied. Then, if a random selec¬ 
tion is done of one of the Sl(Ei) microstates with energy 
Ei and a new microstate is generated by applying a dis¬ 
placement Ar to a random particle, then the probability 
that the system has a new energy Ej is given by 


Similarly, 


P{Ei -A Ej) 


Vij 

NSi{E t ) ’ 


Vji 

NQ(Ej)' 


(3) 

(4) 


Both probabilities can be obtained calculating the rate 
of attempts T i:j to go from level Ei to level Ej. In order 
to do this, two variables are required: 

• The number of times that the system spend in level 
Ei, denoted by z*. 


• The number of times that the system attempts to 
go from level Ei to level Ej, denoted by Zij. 

Variables Zi and z^ can be evaluated according to the 
following steps: 


1. With Ei as the initial state, a particle at random 
is chosen and Zi is then redefined as Zi + 1. 

2. The new energy Ej is evaluated considering that a 
random displacement is applied to the previously 
chosen particle. 

3. If Ej is an allowed energy level, then zy is up¬ 
dated as + 1, independently of the particle’s dis¬ 
placement being accepted or not. The only possible 
restriction is the value of the fixed energy chosen 
to work with, i.e. all cases where Ej < E m i n or 
Ej > E max are discarded. 

4. The probability of an accepted displacement is 1 if 
Tij < Tji , otherwise is Tji/Tij. 

The last condition assures that all levels are able to be 
visited with equal probability, independently of their de¬ 
generacy. The initial values for z, and Zjj can be any 
positive number and after a large number of particle’s 
displacement attempts it is observed that 


and the required ratios are given by 

— tt(Ej) 

Tji fl(Ei) 


(5) 

( 6 ) 


This algorithm is highly efficient to obtain the ratios 
n(Ei)/Q(Ej) (or the entropy differences Si — Sj, accord¬ 
ing to Boltzmann’s relation S(E) = fcIn [f2(£7)]), since 
the number of times that the random number generator 
is used is smaller than those required in MC-NVT sim¬ 
ulations. The efficiency of the method increases if the 
number of allowed levels ( E max — E m i n )/e is decreased. 

Since the inverse temperature /3(E) = 1 /kT is obtained 
from the entropy by deriving it with respect to the energy, 


p (E) = dS/dE (7) 


then it is convenient to express the entropy as a series 
expansion in /3{Ei ), 


S(Ej) — S(Ei) + erj(3(Ei) + ..., 


P{Ej -A Ei) 


( 8 ) 
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where 77 is an integer such that Ej = Ei + rje. The rest 
of the terms in the expansion can be discarded as long as 
N is large enough, obtaining 


1 n (Tij/Tji) 


er] dS 
~k dE 


(9) 


This equation can be used to obtain the inverse temper¬ 
ature as a function of the internal energy. 

In the Ising model the energy changes are well defined 
since they depend on the number of nearest neighbors in 
each lattice site and this number is fixed, and the HPM 
algorithm can be applied straightforwardly. In the next 
section we will explain how this method can also be valid 
for the case of a SW fluid. 


III. THE SW FLUID CASE 

We consider SW systems with attractive ranges A= 
1.1, 1.3 and 1.5, that are typical values required to de¬ 
scribe molecular and complex fluids. Previous to any 
calculation we must establish the allowed values of 77 for 
the SW potential. In the original work by Hiiller and 
Pleimling this is not necessary since for the Ising model 
the energy changes are well defined and depend on the 
number of nearest neighbors, which is a fixed quantity. 
For the SW fluid case, however, different jumps of en¬ 
ergy (i.e., different 77 values ) will occur, with different 
frequencies among them. Nevertheless, the values of 77 
can be chosen by considering the more frequent energy 
jumps performed by the system. This can be achieved 
using a simple NVT Monte Carlo simulation for several 
values of density and A in the limit of infinity temper¬ 
ature, a limit in which all energy changes are accepted. 
Then histograms can be obtained considering the relative 
frequency of the energy changes, as presented in Figure[T| 
for three different densities of the SW fluid with A = 1.5. 
Applying this method for the densities and SW ranges 
values used in this work, accurate results can be obtained 
for the inverse temperature /3(E) when 77 = ±1, ± 2 and 
±3. If we restrict the analysis with 77 = ±1 it is possible 
to obtain reliable values of j3{E) but at the cost of in¬ 
creasing the computing time of the simulations. In any 
case, it is necessary to evaluate proper values for 77 each 
time that a SW system is simulated using HPM. 

Once that 77 values have been obtained, the inverse 
temperatures are given by 

1 3 1 

P ( E t ) = -- ^ - ln(T M+n /T i+7?ii ), 77 ^ 0. (10) 

° 77=-3 V 

Furthermore, from the curves /3(E) = S'(E) the isochoric 
heat capacity c(E) can also be obtained performing a 
second derivative, 



A U/e 


FIG. 1: Histograms for the energy displacements 77 = A U/e 
for the SW fluid with attractive range A = 1.5 and densities 
p* = 0.1, 0.4 and 0.7, from top to bottom. 


In order to illustrate this method, computer simula¬ 
tions for the SW fluid were performed using a unitary box 
with periodic boundary conditions, considering N = 512 
particles and reduced densities p* = pa 3 between 0.1 
and 0.8. In all the cases, the reduced energy u* = E/Ne 
values are restricted between u* nin and u* max i n the su¬ 
percritical region. The number of particle’s displacement 
attempts considered were from N x 10 ^, for the smaller 
energy intervals, to 2.5 x N x 10 7 for the higher energy 
intervals and performing 8 different independent runs. 

The inverse temperature obtained values /?* = e/kT 
were fitted to a second-order degree polynomial on u* 

/3* = ao + aiu* + a 2 u* 2 , ( 12 ) 

The coefficients obtained from these fitted expressions 
are given in tables I-III for the values of density and 
attractive ranges used in this work, as well as the range of 
validity for the fitted expression. In Figure [2] we present 
results for p* = 0.4 and three different SW systems. 


TABLE I: Second-order polynomial fit coefficients for A = 1.1. 


* 

p 

do 

an 

a 2 

* 

“'min 

7/* 

'-'"max 

0.10 

-0.96087 

-14.54988 

-23.13050 

-0.225 

-0.008 

0.20 

-1.17362 

-7.82579 

-6.07659 

-0.410 

- 0.200 

0.30 

-1.35121 

-5.20500 

-2.33237 

-0.700 

-0.320 

0.40 

-1.62688 

-4.10229 

-1.27205 

-0.910 

-0.495 

0.50 

-1.88153 

-3.25117 

-0.65977 

- 1.200 

-0.700 

0.60 

-2.25890 

-2.78357 

-0.39227 

-1.520 

-0.980 

0.70 

-2.86484 

-2.59557 

-0.27284 

-1.880 

-1.380 

0.80 

-3.27924 

-2.10602 

-0.09942 

-2.350 

-1.750 


c(E) 


\my 

S"(E) ' 


( 11 ) 
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FIG. 2: (Color on line) Inverse temperature as a function of 
energy u* for density p* = 0.4 for the SW fluid with attractive 
ranges A = 1.1, 1.3 and 1.5. Black dots denote the computer 
simulation data obtained in this work, and the red solid line 
is a second-order polynomial fit. 



FIG. 3: (Color on line) Excess energies of SW fluids with at¬ 
tractive ranges A = 1.1, 1.3 and 1.5, for the isotherm T* = 2.0 
and several densities p*. Black circles correspond to results 
obtained with the polynomial expressions with coefficients 
given in tables I-III, and red triangles are NVT MC simu¬ 
lation data generated in this work. 


TABLE II: Second-order polynomial fit coefficients for A = 
1.3. 


p* 

a o 

ai 

a 2 

^min 

'U'max 

0.10 

-1.40800 

-6.26420 

-4.01083 

-0.550 

-0.300 

0.20 

-2.01277 

-4.20715 

-1.44496 

-1.000 

-0.650 

0.30 

-2.89267 

-3.80293 

-0.90423 

-1.430 

-1.050 

0.40 

-4.26284 

-3.99128 

-0.73094 

-1.920 

-1.500 

0.50 

-6.41629 

-4.57759 

-0.67761 

-2.450 

-2.000 

0.60 

-9.80357 

-5.57877 

-0.68950 

-3.300 

-2.630 

0.70 

-8.45564 

-3.31269 

-0.21241 

-3.630 

-3.280 

0.80 

-6.62135 

-1.48119 

0.05530 

-4.300 

-3.950 


From the fitted expressions it is possible then to evalu¬ 
ate the energies and heat capacities, for a given temper¬ 
ature T* in the range of validity of the polynomial (fl2l) . 
In Figure [3] we present the energy values for the isotherm 
T* = 2.0; results are compared with conventional MC- 
NVT simulated values, obtained using 864 particles, with 
2.5 x 10 5 cycles required for equilibration and 5.0 x 10 5 
cycles to obtain averaged quantities. 

In Figure [4] results are presented for the reduced iso- 
choric heat capacity, c* = c(E)/Ne and are compared 
with MC-NVT values reported by Largo et al. [25j, ob¬ 
taining a remarkable compatibility between results. 



FIG. 4: (Color on line) Isochoric heat capacities of SW flu¬ 
ids with attractive ranges A = 1.1, 1.3 and 1.5, for several 
temperatures T* and densities p*: A = 1.1 with T* = 1.0, 
1.5 and 2.0 from top to bottom, A = 1.3 with T* = 1.5, 2.0 
and 2.5 from top to bottom, and A = 1.5 with T * = 1.5, 2.0 
and 2.5 from top to bottom. In all cases black circles cor¬ 
respond to results obtained with the polynomial expressions 
whose coefficients are given in tables I-III and red triangles 
are MC-NVT results from reference [25| . 
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TABLE III: Second-order polynomial fit coefficients for A = 
1 . 5 . 


* 

p 

a o 

a i 

a 2 

77* • 
^min 

max 

0.10 

- 1.65393 

- 3.97199 

- 1.56668 

- 1.100 

- 0.574 

0.20 

- 2.94547 

- 3.55764 

- 0.84737 

- 1.795 

- 1.205 

0.30 

- 5.15246 

- 4.16089 

- 0.72716 

- 2.495 

- 1.900 

0.40 

- 8.96461 

- 5.41218 

- 0.74406 

- 3.145 

- 2.654 

0.50 

- 11.88410 

- 5.45835 

- 0.56681 

- 3.825 

- 3.400 

0.60 

- 9.17055 

- 2.83072 

- 0.14575 

- 4.555 

- 4.200 

0.70 

- 7.50855 

- 1.45483 

0.01771 

- 5.295 

- 4.950 

0.80 

- 1.37322 

1.14105 

0.24968 

- 5.975 

- 5.650 


IV. CONCLUSIONS 

In this work we have adapted a microcanonical algo¬ 
rithm, originally developed for the Ising model, to simu¬ 
late thermodynamic properties of fluids whose molecules 
interact via a pair SW potential of variable range. The 
algorithm has the advantage with respect to a MC-NVT 
simulation that a continuous range of temperatures is ob¬ 
tained for a given density, and it is also possible to predict 
accurate results for the internal energy and the isochoric 
heat capacity. Since the method is based on the Ising 


model using discrete values, in the case of SW fluids can 
only be applied to evaluate entropy derivatives with re¬ 
spect to discrete quantities, like the number of particles, 
that would be useful in order to obtain chemical poten¬ 
tials. Since the SW fluid is the key ingredient to study 
other DP systems, we are currently studying its extension 
to these generalized models of fluids. Finally, to imple¬ 
ment a protocol for continuous systems could be possible, 
considering that the probabilities to perform and revert 
each allowed movement in the states space of the system 
are the same, and defining a probability to obtain an 
energy change starting from a given configuration [l 4| . 
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